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ABSTRACT 

The evolution of FRI jets has been long studied in the framework of the FRI-FRII dichotomy. 
The present paradigm consists of the expansion of overpressured jets in the ambient medium 
and the generation of standing recollimation shocks, followed by mass entrainment from 
the external medium that decelerates the jets to subsonic speeds. In this paper, we test the 
present theoretical and observational models via a relativistic numerical simulation of the 
jet s in the radio ga l axy 3C 3 1 . We use the parameters derived from the modelling presented 
by lLaing & Bridle! (2002a b) as input parameters for the simulation of the evolution of the 
source, thus assuming that they have not varied over the lifetime of the source. We simulate 
about 10 % of the total lifetime of the jets in 3C 31. Realistic density and pressure gradients 
for the atmosphere are used. The simulation includes an equation of state for a two-component 
relativistic gas that allows a separate treatment of leptonic and baryonic matter. We compare 
our results with the modelling of the observational data of the source. Our results show that 
the bow shock evolves self-similarly at a quasi-constant speed, with slight deceleration by the 
end of the simulation, in agreement with recent X-ray observations that show the presence 
of bow shocks in FRI sources. The jet expands until it becomes underpressured with respect 
to the ambient medium, and then recollimates. Subsequent oscillations around pressure equi- 
librium and generation of standing shocks lead to the mass loading and disruption of the jet 
flow. We derive an estimate for the minimum age of the source of t > 1. 10 8 yrs, which may 
imply continuous activity of 3C 3 1 since the triggering of its activity. The simulation shows 
that weak CSS sources may be the young counterparts of FRIs. We conclude that the observed 
properties of the jets in 3C 31 are basically recovered by the standing shock scenario. 
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- radio continuum: galaxies 



1 INTRODUCTION 

The morphology and evol ution of jets in low - power radiogalaxies, 
i.e., FRI radiogalaxies in iFanaroff & Rilevsl J 19741) classification, 
have been addressed by a number of auth or s, from theoretical (e.g ., 
iBicknellll 19841 1 994lKomissarovl | 1 990 M 1 1 994 J5e" Youndl 1 9931) . 
numerical (e.g ..|Pe Yound 1 1 9861; iBowman et alJI 1996b a nd obser- 
vational (e.g., IParma et al.l l 20021 : lLaing & Bridle! l2002lbT) points 
of view. FRI jets show bright regions close to the core dimming 
into diffuse emission at kiloparsec scales. The jets in FRI so urces 
are relativistic at the parsec scales (see, e.g.. lLaing|l993L[l996t) . but 
dec elerate betwe en the inner regions and the kiloparsec scales (see, 
e.g.. lLaing|| 1996b . 

The theoretical paradigm for the evo lution of jets in FRI 
sources (Bicknell 1984 iLaingl Il993lll996b comprises the expan- 
sion of an overpressured jet, followed by the generation of shocks 
due to subsequent compressions and expansions around pressure 
equilibrium; the jet is decelerated in these shocks and entraines 
the external medium through turbulent mixing. The last stage of 



the evolution, when the jet has already been decelerated to sub- 
sonic speeds, is dominated by turbulence. Two main processes for 
the entrainment of the ambient medium in the jets have been in- 
voked and studied in the literature: 1) entrainment through mixing 
in a turbulent shear layer between th e jet and the medium , and 2) 
entrainment from stellar mass losses. [Komissarovl dl990 a b) devel- 
oped a theoretical model for the mass entrainment at the jet bound- 
ary and made steady-state numerical calculations, succeeding in 
explaining the main featur es of FRI jets. A series of autho r s have 
concentrated on this topic. |Pe Yound dl986h and iBicknelll d 1994b 
studied the entrainment in jets. |Pe Yound dl993b studied the im- 
portance of a suffici ently dense enviro nment in decelerating jets 
through entrainment. iKomissarovl d 1994b studied the mass load of 
a leptonic jet by stellar winds, concluding that this e ntrainment can 
be ver y important for the process of deceleration. iBowman et al.l 
( 1996) performed a series of steady-state numerical simulations us- 
ing an equation of state for a two-component relativistic gas (see 
also lBowmanll994) and including a term for the entrainment of the 
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stellar mass loss in the code according to the model of Kom issarovl 



2 SETUP OF THE NUMERICAL SIMULATION 



T — I c t> jiJ JinmJITJ, rnm ur <- j 2.1 The model of Laing & Bridle 

Laing & Bridle (2002a b) -LB02a,b from now on- presented & 



a model which accurately describes the kinematics and dynamics 
of the jets in the FRI radiogalaxy 3C 31, mapping the emission 
and magnetic fields of the jets. In LB02a, the observed brightness 
and polarization distributions were fitted by modelling the veloc- 
ity, synchrotron emissivity and ordering of the magnetic field. In 
LB02b, a dynamical model was presented based on the results of 
LB02a and estimates for external pressure and density profiles from 
Chandra ( lHardcastle et alj|2002l) . applying conservation of parti- 
cles, momentum and energy, and assuming that the jets are in pres- 
sure equilib rium with the external medium at large distances from 
the nucleus. lLaing & Bridlelfe004l) studied the validity of so-called 
adiabatic models in explaining the structure of the magnetic and ve- 
locity fields, and the brightness and polarization distributions in the 
jets of 3C 31. In these models, the radiating particles are supposed 
to be accelerated before entering the region of interest and then lose 
energy only by adiabatic processes. The authors concluded that the 
adiabatic models give a good description of the outer regions of the 
jet, whereas closer to the nucleus, the jet shows a non-adiabatic be- 
havior. Their fit to the data is inferior to that of the free models of 
LB02a, but is obtained with fewer free parameters. 

In this paper, we present the results from a relativistic, 
purely hydrodynamical simulation with input parameters taken 
from LB02a,b. It is not clear if the magnetic field in the jets of 
3C 31 is dynamically important, although it seems, from the re- 
sults obtained by LB02a,b, that it is close to equipartition. How- 
ever, we will assume in this work that the magnetic field is not 
important for the dynamics of the problem. This assumption re- 
mains to be checked by future relativistic magnetohydrodynamic 
simulations. Our aim is to compare the results from the simula- 
tion with those from observations and modelling. We also address 
the problem of the jet evolution and young counterparts to FRI 
jets. The nature of compact radio-sources and their relation with 
large scale jets has been studied during the last years. Compact 
Symmetric Objects (CSO), as a morphological subclass of Giga- 
hertz Peaked Spectrum sources (GPS) and Compact Steep Spec- 
trum sources ( CSS), are understood as the young coun terparts of 
FRII jets (e.g jFanti et al.ll995l ; |Perucho &~M arti 2002^). However, 
there is still so me debate on the com pac t radio sources that could be 
the yo ung FRIs. lDrakeetal.ld2004ri and lKunert -Bairaszewska et al.l 
d2005h have observed low-pow er CSS from the IR AS sample and 
the FIRST survey, respectively. iDrake et"al] ^20041) proposed these 
weak CSS sources as the possible young counterparts of FRI jets. 
We study this possibility in this work. 

To our knowledge this is the first work to present 
long term relativistic simulations of FRI jets. In contrast, the 
evolution of F RII jets has been addres s ed by several au- 
thors (see, e.g.. iKomissarov & Faliel Il998l : IScheck et all |2002| ; 
iKrause & Camenzindll2003l : lKrausell2005l) . Our work covers a pa - 
rameter space complementary to that used in IScheck et alj J2002h . 
who discussed the dependence of the morphology and dynamics of 
jets on their composition, using the same equation of state as used 
here. 

The plan of the paper is as follows. In Section [2] we present 
the setup of our simulation. In Section [3] we summarize the main 
results of the simulation, and we compare these results with new 
simulations performed in order to test different evolutionary sce- 
narios by changing the initial conditions. In Section [4] we discuss 
the main results of the work, and we present the conclusions of the 
paper in Section[5] 



The numerical simulation presented in this paper is based on the 
work by Laing & Bridle (LB02a,b), as stated in the Introduction. 
In this paragraph we review their main assumptions and conclu- 
sions. In the modelling presented in LB02a the jet and counter-jet 
are assumed to be identical, antiparallel, axisymmetric and station- 
ary, and the differences between them are assumed to be due to 
relativistic aberration. The models are designed to fit the observed 
brightness and polarization distributions (assuming optically thin 
emission), taking into account Doppler boosting effects and relying 
on simple prescriptions for the variations of the flow velocity, syn- 
chrotron emissivity an magnetic-field structure within the jet. The 
rotation of the line of sight relative to the magnetic field structure 
between the emitted and the observed frames is also taken into ac- 
count. The model focuses on the region enclosing the inner 12 kpc 
of the jets, which is split in three parts: the inner region (from to 
1.1 kpc, or from to 2.5 arcsec), the flaring region (from 1. 1 to 3.5 
kpc, or from 2.5 to 8.3 arcsec) and the outer region (from 3.5 to 12 
kpc, or from 8.3 to 28.3 arcsec). The linear distances are calculated 
assuming a Hubble constant Ho = 70kms _1 Mpc -1 , taking the 
redshift of the parent galaxy of 3C 31 (NGC 383, z = 0.0169) 
and a viewing angle of 52°. This modelling has b een extended 
to other FRI sources ( e.g., NGC 315 and 3C 296. ICanvin et al.l 
l2005l ; lLaing et alj|2006h . In LB02b, the authors present a dynami- 
cal model for the jet, based on: 1) the results obtained in LB02a, 
2) estim ates of pressure and density profiles from Chandra and 
ROSAT ( lHardcastle etafll2002l : iKomossa & BoehringeJll999h . 3) 
the conservation of particles, momentum and energy, 4) the as- 
sumption that the jets are in pressure equilibrium with the external 
medium at large distances from the nucleus, and 5) the momentum 
flux being n = $/c, where $ is the energy flux (a good approx- 
imation for light, hot and/or relativistic jets). The model is quasi- 
one-dimensional, as, although the widening of the jet is considered, 
only the axial velocities are used in the analysis. The authors con- 
clude that after the expansion of the jet in the inner region, the jet 
flow is overpressured and decelerated at the beginning of the flaring 
region. In this region, they find local minima of pressure and den- 
sity and maxima in the Mach number and entrainment rate. At the 
end of the flaring region, the jets are slightly underpressured but 
close to pressure equilibrium with the ambient medium. Changes 
in the outer region are smooth, with almost constant density and 
monotonically increasing entrainment rate. The Mach number in 
the outer region is ~ 1 — 2. Comparison of the jet density with the 
number of radiating particles required by the observed synchrotron 
emissivity lead the authors to conclude that the jet is probably ini- 
tially composed of a pair plasma (e~ — e + ), and mass-loaded with 
baryons inside the galaxy. This mass load is attributed dominantly 
to stellar wind material. The exact origin of the majority of the gas 
entrained by the jets is, however, still uncertain. 



2.2 Setup 

The numerical simulation was performed using a finite-difference 
code based on a high-resolution shock-capturing scheme which 
solves the equations of relativistic hydrodynamics i n two dimen- 
sions written in conservation form jMartf et al]|l997h . The code is 
parallelized using OMP directives and it has been modified in or- 
der to include the equation of state for r elativistic Boltzmann gases 
dSvngdl 19571 : iFalle & Komissarovll 19961 . , and Appendix A in this 
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paper) with the routines used in the simulations of IScheck et al.l 
The use of this equation of state allows us to distinguish 
electrons, positrons and protons in the simulation. Under the as- 
sumption of charge neutrality, only one conservation equation (for 
example, that for the evolution of the leptonic density) must be 
added in order to know the composition of the fluid at each com- 
putational cell. The integration of this extra equation together with 
the use of the Synge equation of state (involving the computation 
of Bessel functions) increases the computational time per iteration 
by ~ 50 % with r espect to the case o f the one component, ideal gas 
equation of state JScheck et afll2002l) . 

The equations of relativistic hydrodynamics in conservation 
form, which are solved by the code, assuming axisymmetry and in 
two-dimensional cylindrical coordinates (R,z), are the following 
(we use units in which the speed of light is set to 1): 



dV 1 dRF R d¥ z _ s 

dt + R OR + dz ~ " 

with the vector of unknowns 



V={D,Di,S r ,S z ,t) t 



(1) 



(2) 



fluxes 

F R = (Dv R ,D l v R ,S R v R +p,S z v R ,S R -Dv R ) T , (3) 

(4) 

and source terms 

(5) 



F z = (Dv z ,Div z ,S R v z ,S z v z + p,S z -Dv z ) T 



S = (0,0,p/R + g R ,g z ,v R g R + v z g z ) T 



The five unknowns D, Di, S R , S z and r refer to the densi- 
ties of five conserved quantities, namely the total and leptonic rest 
masses, the radial and axial components of the momentum, and the 
energy (excluding rest mass energy). They are all measured in the 
laboratory frame, and are related to the quantities in the local rest 
frame of the fluid (primitive variables) according to 



D — pW , 

D t = piW , 

S R ' Z = phW 2 v R ' z , 

r = phW 2 - p - D 



(6) 
(7) 
(8) 
(9) 



where p and pi are the total and the leptonic rest mass density, 
respectively, v R ' z are the components of the velocity of the fluid, 
W is the Lorentz factor (W = 1/ \/l — v i v i , where summation 
over repeated indices is implied), and h is the specific enthalpy 
defined as 



h = 1 + e +p/p 



(10) 



where e is the specific internal energy and p is the pressure. Finally, 
g , g z appearing in the definition of the source vector, eq. lO, are 
the components of an external force that keeps the atmosphere in 
equilibrium (see below). The system (TJl is closed by means of the 
Synge equation of state described in Appendix A. 

The code also integrates an equation for the jet mass fraction, 
/. This quantity, set to 1 for the injected beam material and oth- 
erwise, is used as a tracer of the jet material through the grid and 
allows to study phenomena like the entrainment of ambient mate- 
rial in the beam and the mixing in the cocoon. 

The medium in which the jet is injected consists of a decreas- 



ing density atmosphere composed of HydrogerQ. The dynamical 
equilibrium of the atmosphere is attained by introducing an exter- 
nal force which compensates initial pressure gradients in the radial 
and axial d irections. The profile f or the number density of such a 
medium is dHardcastle et alj|2002h : 

c/2 / 2 \ -3/8o«m,g/2 



t(r) 



1 + - 



"3,3a 



,(H) 



where r is the spherical radial coordinate, n c and n g are the core 
densities of the galaxy and the surrounding group, r c and r g are 
the radii of those cores, and P a tm,c and /3 a tm,g are the exponents 
giving the profile for each medium. The temperature profile is: 



T = T c + {T g -T c )— {r<r m ) 



T — T a 



(r ^ r m ), 



(12) 



with r m = 7.8 kpc. The pressure is derived from the following 
equation of state: 



k B T 

1 

pX 



(13) 



where p is the mass per particle in a.m.u. (p — 0.5 in our case, cf. 
0.6 in LB02b), X is the abundance of hydrogen by mass (X — 1 
here, cf . 0.74 in LB 02b). In TableQjwe reproduce Table 1 in LB02b, 
where the parameters for the equations above are listed. The initial 
profiles for ambient medium pressure and density along the axis of 
the jet are plotted in Fig.Q] 

The numerical code units are the jet radius (Rj), the speed 
of light (c) and the density of the ambient medium at injection 
(Pa,c). Thus, the appropriate unit transformations are performed in 
the code from physical to code units. 

The simulation was performed in cylindrical coordinates with 
axial symmetry, i.e., only one half of the jet is computed. The grid 
involved 2880 x 1800 cells, with a resolution of 8 cells/Rj in the 
axial direction up to 300 Rj, plus an extended grid with geomet- 
rically increasing cell size (up to 450 Rj), and 16 cells/Rj in the 
radial direction up to 100 Rj, plus an extended grid (also with ge- 
ometrically increasing cell size) up to 200 Rj. Outflow boundary 
conditions were used at the end of the grid in the axial direction and 
also far from the jet axis in the radial direction. Reflecting bound- 
ary conditions were taken for the jet axis in order to account for the 
cylindrical symmetry. 

Injection of the jet in the atmosphere is done at r = 500 pc, 
the point where the modelling of the jet in LB02b starts. The ra- 
dius of the jet at that point is calculated from the opening an- 
gle of 6.7° given for the jet in LB02a (Rj = 60 pc). The uni- 
form grid is then 17.5 kpc x 6 kpc, and, with the extended grid, 
26.25 kpc x 12 kpc. 

The jet is injected with a speed Vj = 0.87 c (jet Lorentz fac- 
tor, Wj ~ 2), internal relativistic Mach number Mj = 2.5, tem- 
perature Tj = 4.1 10 9 K, density ratio with the external medium 
rj = 1. 10 -5 , purely leptonic composition (Xi — 1.0), and over- 
pressured by a factor 7.8 with respect to the ambient medium. 
In Table [2] we give the complete list of parameters. The parame- 
ters given for the ambient medium at injection are calculated for 
r — 500 pc (this is the reason for differences between Table Q] 

1 LB02b use the standard composition with 74% hydrogen, but this treat- 
ment would require the inclusion of new populations of particles (in order 
to account for the remaining 26%) in the code, involving longer computa- 
tional time, so that we discarded this option. 



4 M. Perucho & J. -MA Marti 



Ambient density Ambient pressure 




5.0x10 3 1.0x10 4 1.5x10 4 2.0x10 4 5.0x10 3 1.0x10 4 1.5x10 4 2.0x10 4 

r(pc) r(pc) 

Figure 1. Initial profiles of ambient rest mass density (a) and pressure (b). 
Table 1. Ambient medium parameters (Table 1 in LB02b). Subscripts c and g refer to the galaxy core and to the surrounding group, respectively. 



Component 


Central density 


Form factor 


Core radius 


Temperature 


Galaxy 


n c = 1.810 5 m- 3 


Patm,c = 0.73 


r c = 1.2 kpc 


T c = 4.910 6 fs: 


Group 


n g = 1.910 3 m- 3 


Patm,g = 0.38 


r g = 52 kpc 


T g = 1.7 W 7 K 



Table 2. Table of parameters used in the simulation. The one-dimensional the authors make a stud y of the influence of jet composition on its 

velocity estimation given in the Table stands for the theoretical advance ve- long term evolution, is equivalent to the jet energy flux defined in 

locity of the jet, computed using the equation derived in Marti et al. (1997) L B02b. Comparing the value of Lk in in our simulation to that used 

in lScheck et al.l d2002t) we use L^ in ~ 10 37 W, on the upper end 
of the power distribution of FRI sources, whereas in their work the 
jets have Lkin ~ 10 39 W, appropriate for an FRII source. 

3 RESULTS 
3.1 Evolution 

3.1.1 Numerical Results 

The simulation was run for ~ 2000 hours (~ 83 days) on eight 
processors in the SGI Altix computer CERCA, at the Universitat 
de Valencia, corresponding to a source lifetime of 7.26 10 6 yrs. At 
the moment when it was stopped, the bow shock of the jet was 
located at ~ 14.5 kpc from the injection point in the grid, i.e., at 
~ 15 kpc from the source. The large amount of computational time 
needed is a consequence of the small advance speed of a FRI-like 
source and of the numerical effort invested in the computation of 
the relativistic equation of state. 

Fig. [2] shows panels of the logarithm of rest mass density and 
Lorentz factor at different times of the first stages of evolution. The 
compact phase, defined as that in which the source has a linear size 
smaller than 5 kpc can be divided into two main epochs, based on 
these maps: 1) the CSO-like phase, when the source is smaller than 
1 kpc, and 2) the weak CSS-like phase, when the source is between 
1 and 5 kpc long. During the CSO-like phase (top panels in Fig. [2)1, 
and Table |2). The numbers in this table give an energy flux for the jet shows a large opening angle, due to overpressure with re- 
the jet of $ ~ 10 37 W, which is very close to the value given in spect to the ambient medium and a strong Mach disk at its head. 
LB02b (<!> = 1.1 10 37 W), the difference being due to approxima- Hence the morphology of the source is dominated by a short and 
tions in variables considered. The quant ity Lkin (kin e tic l uminosity featureless beam and a strong hot-spot downstream of the termi- 
of the jet) used for the simulations in IScheck et al] d2002h . where nal Mach shock. During the weak CSS-like phase (mid and bottom 



for a pressure-matched jet propagating in one dimension (i.e., without side- 
ways expansion) through an homogeneous medium. 



Velocity (vj) 


0.87 c 


Mach number (Mj ) 


2.5 


Temperature (T,- , jet) 


4.1 10 9 K 


Temperature (T c , ambient 1 ) 


5.710 6 K 


Temperature (T g , ambient 2 ) 


1.710 7 K 


Density (pj, jet) 


3 10- 27 kg/m 3 


Density (p a ,c, ambient 1 ) 


310- 22 kg/m 3 


Density ratio (rf) 


10" 5 


Leptonic number (X; , jet) 


1.0 


Specific int. energy (sj, jet) 


1.54 c 2 


Specific int. energy (e a . c , ambient 1 ) 


1.5710" 6 c 2 


Specific int. energy (e a , g , ambient 2 ) 


4.6910" 6 c 2 


Pressure (Pj , jet) 


6.9110- 6 p a ,cc 2 


Pressure (P a , c , ambient 1 ) 


8.8410- 7 p a , c c 2 


Pressure (P a .g, ambient 2 ) 


3.0710- 8 p a ,cc 2 


Pressure ratio (Pj / P a ,c) 


7.8 


Adiabatic exponent (Vj , jet) 


1.38 


Adiabatic exponent (r a , ambient) 


1.66 


ID velocity estimation (ui ) 


9.910~ 3 c 


Time unit (Rj /c) 


60 pc/c ~ 195 yrs 



1 and 2 stand for values in the ambient medium at the injection and point 
furthest from injection in the grid, respectively. 
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panels in Fig.[2]l, the Mach disk is weaker, as the jet is slowed down 
in recollimation shocks. The whole structure of the jet appears to 
be much more irregular during the latter phase. The transition be- 
tween the phases occurs in the course of the evolution, when the 
terminal Mach shock disappears for the first time leading to a con- 
ical shock, at t — 4.3 10 yrs, when the linear size of the jet is 
in the range 1 — 1.5 kpc, i.e. in the transition between CSO and 
CSS phases of young radio sources. The flow dynamics behind the 
shock then changes affecting t he global evolution o f the jet. This 
behaviour was already noted bv lScheck etai] j2002l) in the context 
of the early evolution of FRII jets. It is important to remark that 
the axial symmetry imposed on the simulation affects the internal 
structure of the jet (internal and terminal shocks) and hence the 
jet dyn amics and advance speed as discussed, e.g., by lAlov et al.l 
However, for the short times/sizes involved in the previous 
discussion, the hypothetical three dimensional effects are expected 
to be negligible. 

Fig. [3] shows the evolution of several quantities versus time 
during the simulation. Fig.[3ji shows the position of the bow shock 
as a function of time. The curve is consistent with a constant- 
velocity expansion up to t — 4.5 10 6 yrs (corresponding source 
size in the axial direction, 9.5 kpc), and with a decelerating expan- 
sion afterwards. Despite the variations, the advance velocity of the 
bow shock, in Fig.[3j:, shows this trend. Initially, the advance speed 
is Vb s ~ 7 10 _3 c, close to the one-dimensional estimate of the head 
of the jet for homogeneous ambient medium, vjf = 9.9 10 -3 c 



(Table [2), but by the end of the simulation, it has decreased to 
vts < 6 10~ 3 c. The duration of the phase with constant advance 
speed is longer than in the case of simu lations with a uniform am- 
bient medium (e.g.. IScheck et al]|2002h because the deceleration 
caused by the widening of the jet along the evolution is now ap- 
proximately balanced by the decrease of the inertia of the ambient 
medium with distance. As stated in the previous paragraph, three 
dimensional effects can modify the change of jet advance speed 
with time. In particular, the one dimensional phase can be short- 
ened and phases of acceleration can appear occasionally when the 
terminal shock changes from planar to oblique. Hence the results 
discussed in this paragraph must be treated with caution. The pres- 
sure behind the bow shock, in Fig. |3j), drops and oscillates until 
the end of the simulation following the behaviour of the ambient 
pressure although with a larger instantaneous slope. The effect of 
deceleration of the bow shock and the increase of temperature in 
the ambient medium at larger distances to the source combine to 
give a mild reduction of the Mach number (Fig.[3jl) along the evo- 
lution, although the shock is still superson ic (Mts ~ 2 — 3 ) by 
the end of the simulation. Recent works bv lKraft et all (2003) and 
ICroston et ail {2007) have detected the existence of shock waves 
surrounding the radio lobes of the relatively young (t ~ 2 10 6 yrs) 
FRI jets of Centaurus A and NGC 3801, respectively. The shock 
waves are revealed by shells of hot interstellar gas emitting in the 
X-rays. The authors derive Mach numbers between 3 and 8 for 
these shocks. The ages and Mach numbers obtained in our simu- 
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Figure 3. Evolution with time of the position (a), pressure (b), instantaneous velocity (c) and Mach number (d) of the bow shock in the axial direction. The 
bow shock position is selected as the first point of the grid (starting from its end) where the speed is larger than 10~ 4 c. The dashed line in panel a indicates 
the best fit to the first half of the evolution (z oc 0.0022 t[yrs] pc); deceleration at the last stages is clear. The bow shock pressure is taken as the value in 
the first maximum, also searched from the end of the grid. The dotted line in panel b shows the pressure of the ambient medium in the region where the bow 
shock is located at the given time. The instantaneous velocity of the bow shock is computed with the discretized derivative of position with respect to time at 
each instant, and Mach number of the bow shock propagating in the ambient medium is computed with that advance speed and the mean value of the sound 
speed in 16 cells (2 Rj = 120 pc) ahead of the bow shock position. Bow shock pressure, velocity and Mach number curves have been smoothed using an IDL 
routine (with a 10 point smoothing) in order to avoid the noise coming from the numerical derivative of the bow shock position. 



lation are in agreement with the results reported in those papers. 
We will further discuss this issue in the next Section. 

The propagation of the shock leaves behind a region of 
shocked ambient material through which the jet propagates. In the 
case of powerful jets, the shock is so strong and its expansion so 
fast that the region encompassed by the bow shock is almost evacu- 
ated with most of the matter concentrated in a thin shell of shocked 
ambient material behind the shock. The evacuated region is con- 
tinuously fed by the jet to form an extended cocoon. In the case 
of weak jets like the one considered here, the bow shocks are cor- 
respondingly weaker and the dense region behind the shock, cor- 
respondingly wider. Hence the cavity/shell division suitable for 
the shocked regions surrounding powerful jets transforms into a 
cocoon/shocked-ambient-medium division. This structure is easily 
seen in the color maps of rest-mass density and tracer at the end of 
the simulation, to be discussed in the next section. 

Fig. [4] shows the evolution of the jet length (Fig. [4^) and the 
jet radius (Fig.|4p) with time. Both quantities are calculated for dif- 
ferent values of the jet mass fraction. The spreading of the lines in 



the jet length plot for t > 4.5 10 yrs results from the entrainment 
of ambient material in the head of the jet which causes the decel- 
eration in the jet advance mentioned in the preceding paragraphs. 
The values of the jet radius for the different mass fractions diverge 
right from the start of the simulation, implying progressive jet strat- 
ification, shearing and mixing, as the jet propagates outwards. The 
aspect ratio of the jet increases for larger values of the jet mass 
fraction and is always larger than the aspect ratio of the shocked 
region. 



3.1.2 Comparison with Analytical Models 

iBegelman & Cioffil d 19891) developed a simple model for the evolu- 
tion of the cavities/cocoons surrounding powerful jets in homoge- 
neous ambient media under the assumptions that the jet propaga- 
tion speed and the power injected into the cavity are independent 
of time, and that the pressure of the external medium is negligi- 
ble. The model can still be applied to describe the evolution of the 
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shocked regions surrounding weak jets. Then, the mean pressure in 
the shocked region is given by 

Ps « (14) 

where L s , Vbs and A s are the power injected into the shocked re- 
gion from the jet through the terminal shock, the bow shock speed 
in the direction of propagation of the jet, and the transversal cross- 
section of the shocked region (i.e., irR 2 , R s being the radius), re- 
spectively. The pressure in the shocked region causes it to expand 
with velocity R s according to 

Ps = PaRl, (15) 

where p a is the ambient density. This implies 1/R S oc R s , and 
hence 

R s oct 1/2 , P s oci _1 , l s /R s <xt 1/2 (16) 

(where in this last expression, l s stands for the longitudinal size of 
the sho cked region). 

In lScheck etai] J2002I) . and also in the context of powerful 
jets, the authors developed a simple extension of Begelman & 
Cioffis's model to account for the secular deceleration of the jet 
advance speed due to the expansion of the jet cross-section. Ac- 
cording to this model, Vbs oc t a , and 

R a oc t 1/2 ~ a/ \ Psoct- 1 -"/ 2 , ls/R s oct 1/2+5a/4 . (17) 

In the simulations presented in that paper, a value of a ~ — 1/3 
was found and, accordingly, 

R s oct 7/12 , PsOct~ 5/6 , L/Rs oct 1/12 . (18) 

In the present work, the model is generalized to consider the 
expansion of the shocked region through an ambient medium with 
decreasing density. A power law fit of the ambient density (see 
Fig. [TJ with respect to the distance to the source gives p a oc RJ 1 . 
With this and taking a ~ —0.1 (from the simulation), we finally 
have 

R s (xt - 7 , PsXt' 1 - 3 , l s /R s (xt - 2 . (19) 

These results are remarkably consistent with those derived from 
the simulation (R s oc t°- s , P s oc t~ 13 , l s /R s oc t 2 ). Fig. [5ji 
displays the evolution of the pressure in the shocked region versus 



time, together with the t" 1 (as in Begelman & Cioffi and Scheck 
et al. models) and i" 1 ' 3 fits. At this point, we must remark that in 
the simulation we used an open boundary condition in the symme- 
try plane at the jet basis that allows the leakage of gas and, con- 
sequently, a faster pressure drop in the shocked region. A crude a 
posteriori analysis does not allow us to rule out the possibility that 
the gas internal energy that left the shocked region through the open 
boundary is of the same order (or ~ 1/10) than that remaining in 
it. However, the agreement between the results from the model dis- 
cussed above and those derived from the simulation suggests that 
the flow through the open boundary is effectively negligible. De- 
spite the continuous decrease of the pressure, the shocked region is 
still overpressured by a factor of 4 with respect to its environment 
at the end of the simulation. This fact migh t be alleviated w ith the 
introduction o f cooling processes. However, iKraft et al.l d2003l) and 
ICroston et al.l fc007h also derived strong overpressure (as much as 
two orders of magnitude) of the X-ray emitting shells in Cen A and 
NGC 3801 with respect to the surrounding interstellar media. This 
is in agreement with our results that give a pressure in the shocked 
region at t ~ 2 10 6 yrs (the approximate age of the jets in Cen A 
and NGC 3801) more than one order of magnitude larger than that 
of the ambient medium. Finally, the plot of shocked region radius 
versus linear size, in Fig.|5p, shows a self-similar growth for the late 
stages of the simulation (with an aspect ratio for the shocked region 
of « 2.7). The deceleration of the jet advance for t > 4.5 10 6 yrs, 
produces a change in the evolution of the aspect ratio of the shocked 
region that by the end of the simulation has decreased to a value of 
2.6. 

Fig.[6]show the evolution of the mean pressure (Fig.|6ji), den- 
sity (Fig.|6p) and temperature (Fig. [(J) in the cocoon as a function 
of time. Remarkably, the cocoon temperature is almost constant in 
the long term evolution. Also interesting is the fact that the pres- 
sure in the cocoon follows a similar evolution to that in the whole 
shocked region (see Fig.[5ji). This last fact can be explained taking 
into account that the sound speed in the cocoon/shocked-ambient- 
medium region, is about one or two orders of magnitude larger than 
its expansion velocity, hence allowing for an almost instantaneous 
(with respect to the dynamical time scale) adjustment of the pres- 
sure. An example of this is seen in the pressure map at the end of 
the simulation to be discussed in the next section. 

The isothermal evolution of the cocoon can be explained by 
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Figure 5. The panels show the pressure in the shocked region versus time (a) and the radius of the shocked region versus its length (b). The pressure is 
computed as the mean pressure in all the cells with shocked ambient medium. The dotted line shows the evolution as oc t~ 1 and the dashed line gives the best 
fit of the curve t -1 - 3 . The radius is calculated as the mean radius of the bow shock, with this position determined by a lateral motion larger than 10 — 4 c, see 
caption of Fig. [3] The length of the shocked region used in this plot as the abscissa coincides with the bow shock position shown in the panel a of Fig. [5] 
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Figure 6. Mean pressure (a), density (b) and temperature in the cocoon as a function of time (c). The cocoon is defined as the region with tracer values between 
0.01 and 0.3. 



A numerical simulation of the evolution and fate of a FRI jet. The case of 3C 31 9 



assuming that the almost self-similar evolution of the shocked re- 
gion extends also to its internal structure. As discussed in the pre- 
vious paragraph, for the pressure in the cocoon, P c , we have, at any 
time, 

P c ~P 3 oc-^-. (20) 

The density in the cocoon, p c , can be estimated, in a similar way, 
as the quotient of total mass injected in the cocoon, J c t, and its 
volume. If Ac is the transversal cross-section of the cocoon, then 



On the other hand, for the temperature in the cocoon, T c , we have 
T c oc P c / pc, and, from the two previous equations, 

T.oc^- (22) 

JC 

Now, the independency of T c follows if the fluxes of energy and 
mass through the jet terminal shock are constant and the evolution 
of A c and A s with time is the same. 

3.2 Fate 

Figs. 17191 show several maps of different quantities at the end of 
the simulation (t — 7.26 10 6 yrs). The morphological features ob- 
serve d in the panels are those of a typical jet (see, e.g jMarti et al.l 
[199% : a bow-shock propagating through the unperturbed ambient 
gas (clearly seen in the top panel of Fig. [7j displaying the loga- 
rithm of the rest-mass density), a cocoon composed of mixed jet 
and ambient matter, and the jet itself propagating inside. The co- 
coon, formed by mixed jet material, is better observed in the bot- 
tom panel of Fig. [7j displaying the jet mass fraction. Both maps 
show how the jet initially expands up to a distance z ~ 1.5 kpc 
from the sourcq 2 ] and then recollimates and oscillates until it is dis- 
rupted, i.e., until the ambient medium material reaches the jet axis 
(at 2 ~ 4.5 kpc), due to entrainment of the external medium. Fig.[8] 
shows the flow Lorentz factor (top panel) and the axial velocity 
(bottom panel). A relativistic jet (with maximum Lorentz factor of 
5.3) is seen up to the disruption point at z ~ 4.5 kpc. After this 
point, strong deceleration of the flow associated with mass loading 
of the jet is observed in both maps. The axial velocity plot reveals 
a mildly relativistic backflow in the cocoon. However, these high 
sp eeds can be an art ifact of the imposed axisymmetry, as discussed 
in lAlov et alj fl999h on the basis of 3D simulations of relativistic 
jets. The use of open boundaries on the symmetry plane at the jet 
base could also maintain artificially large pressure gradients along 
the cocoon leading to high speed backflows (see the discussion in 
sect. 13X21 Fig.|9]shows maps of the logarithm of temperature (top 
panel) and pressure (bottom panel). The bow shock, cocoon and jet 
structure are also clearly seen in these maps. As can be seen in the 
top panel of Fig. [9] the ambient temperature increases by a factor 
of a few behind the bow shock and by two orders of magnitude in 
the cocoon due to the mixing with hot plasma injected through the 
terminal shock at the head of the jet. Also in this panel, Kelvin- 
Helmholtz instabilities arising in the interface between the outer 
subsonic backflow in the cocoon and the shocked ambient medium 
are observed. Finally, the bottom panel of Fig. [9] shows a shocked 

2 Throughout this section and in the discussion of results, z will refer to 
distances to the galactic source, located 500 pc away from the injection in 
the grid. 



region with an homogeneous pressure distribution and a series of 
waves emanating from the outer part of the cocoon. 

The profiles of several variables along the jet axis (rest-mass 
density, pressure, jet mass fraction, Lorentz factor, axial velocity 
and temperature) at the end of the simulation are plotted in Fig.UOl 
Fig.[TT]shows averaged values of rest-mass density, pressure, Mach 
number and flow Lorentz factor along the jet at the same time. From 
these plots and the maps in Figs. |7|9| we can obtain a general picture 
of the evolution of the jet in terms of its dynamics. Close to the in- 
jection point, the flow can be considered as steady. The jet initially 
expands in the ambient density gradient. The expansion accelerates 
(Fig. [8), rarefies (Fig. [7J and cools (Fig. |9]l the flow. This is seen 
in more detail in the cuts along the jet axis shown in Fig.QJj] After 
a short distance, in which the variables remain constant, a strong 
adiabatic expansion produces a fast decrease in density, pressure 
and temperature on the jet axis, up to z = 1.5 — 2.0 kpc. At the 
same time, there is a strong acceleration seen in the plots of Lorentz 
factor and axial velocity. The jet overpressures again in a stand- 
ing shock at z ~ 2 kpc from the source. The shock is seen in the 
maps and plots of density, pressure and temperature as a sudden 
increase and, in those of velocity and Lorentz factor, as a strong 
deceleration of the flow. The overpressure of the jet with respect 
to the ambient (see the pressure map in Fig. [9] and the plots in 
the Figs. 1 10b and II lb) results in a new phase of expansion. Each 
expansion is followed by a recollimation shock that produces sig- 
nificant deceleration of the jet (see Fig.[8]and Figs. HOB andllOfe). 
Up to three recollimation shocks (at z ~ 2 kpc, z ~ 3 kpc and 
z ~ 4.5 kpc) are observed before the jet is disrupted. These shocks 
can be clearly identified in the Lorentz factor and velocity plots in 
Fig. [10] (panels d and e). The planar/conical shape of these shocks 
could be an artefact of axisymmetry, but at the distances from the 
central source at which these shocks are formed (i.e. a few kpc), 
3D effects are probably small. After the first shock, the jet is less 
overpressured with respect to the ambient medium and it expands 
in an atmosphere with a smoother density gradient. This makes the 
next standing shocks milder (see the pressure plots in Fig.llOb and 
Fig. II lb). At z > 6 kpc (i.e., behind the terminal shock), the jet 
is overpressured with respect to the ambient (Fig. [TTb ) due to the 
heating of the flow in shocks, as seen in Fig.llOf. 

The Lorentz factor is between 2 and 5.3 in the whole section of 
the jet up to the first recollimation shock, at z ~ 2 kpc (see Fig. [8}. 
After this shock, the flow is relativistic in filaments far from the 
axis, as the portions of the jet closer to the axis are decelerated by 
the shock. The averaged Lorentz factor (Fig. II IB ) decreases from 
4.3, after the first adiabatic expansion of the jet, to a value of 1.15 
at the disruption point, which we define as the point in which am- 
bient material reaches the axis (see the jet mass fraction map in 
Fig. [7] the flow Lorentz factor map in Fig. [8] and the correspond- 
ing axial cuts in Figs. 1 10b and 11 OB). The internal Mach number of 
the flow (Fig. II lb ) also decreases along the jet up to the terminal 
shock where the flow becomes subsonic. After the terminal shock, 
the flow is only slightly relativistic, with velocities around 0.5 c. It 
is also remarkable that in the cocoon a backflow with mildly rela- 
tivistic speeds is established right down to the central regions of the 
parent galaxy. 

The jet remains well collimated and relativistic up to the dis- 
ruption point. However, after this point, the subsonic character of 
the flow triggers the proces s of mixing wi th the ambient medium, 
as already pointed out by iBicknelll dl995l) . The rest-mass density 
plots (Figs. 1 10b and II lb ) show the increase of this variable as 
the jet entrains ambient material for z > 4.5 kpc. The position 
of the disruption point oscillates throughout the simulation from 
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Figure 7. Logarithm of rest mass density (top panel) in kg/m s and jet mass fraction (used as a tracer for the jet plasma; bottom panel) at the last frame of 
the simulation in this work. The locations of the first recollimation shock and the terminal (or bow) shock are indicated in the rest mass density map, and the 
disruption point is indicated in the jet mass fraction map. Coordinates are in parsecs. 



z ~ 2.5 — 3.5 kpc to z ~ 6 — 7kpc. This non-monotonic ad- 
vance is due to the complex dynamics of the terminal shock, which 
changes periodically from planar to conical, as already noted in 
previous long-term simulations of jets (see, e.g.. lMartf et all 1997c 
Schec k et al. 2002). In the phases in which the disruption point 
moves backwards, the further injected plasma accumulates in the 
region between injection and disruption, thus increasing the pres- 
sure in the jet, that ultimately bursts, expanding the unmixed jet 
outwards and bringing the disruption point farther from injection. 

Fig. [9] shows that the temperature in the cocoon is of the 
same order, but larger than that of the jet (T ~ 10 9 - 10 10 K, 



cf Tj —4.1 10 9 K at injection). In addition, the gas in this high 
temperature region can be up to an order of magnitude denser than 
the jet (top panel in Fig. [7]l due to the loading with baryons from 
the ambient (jet mass fraction, / < 0.2). However, despite this 
large baryon load, the number of leptons is still between 20 and 
200 times larger than the number of baryons in this region. Leptons 
are heated in shocks inside the jet (see top panel of Fig. [9} and at 
the terminal shock. The high temperatures achieved in the cocoon 
are a consequence of the small amount of baryons. The adiabatic 
exponent in the cocoon remains close to the relativistic limit of 4/3. 

In a recent paper, Kino, Kawakatu and Ito (2007) conclude 
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Figure 8. Lorentz factor (top panel) and axial velocity (bottom panel) in units of c at the last frame of the simulation in this work. Coordinates are in parsecs. 



that the temperature of the cocoon gas for an FRII jet with Lorentz 
factor 10 could be of about 1.2 10 10 K, a value close to that ob- 
tained in our simulation. In the case of FRII jets, the particles are 
heated in the hot spot, the heating depending on the strength of the 
shock. In a electron-positron FRI jet like that simulated here, the 
particles are heated in the strong recollimation shocks along the jet. 
This effect can compensate the lack of heating at the terminal point 
of the FRI jet due to its low power. Thus, the temperature in the 
cocoon remains high as long as the pollution by baryons is low. 



3.3 Additional simulations 

The modelling of 3C 31 presented in LB02b divides the inner 
12 kpc of the jet into three regions: inner, flaring and outer region. 
The authors suggest that the boundary between the inner and flar- 
ing regions (at 1.1 kpc from the source) consists of a discontinuity 
in velocity, density and press ure, which is the cause of the sudden 
increa se of radio emission. In Canvin & Laing ( 2 0041) . ICanvin et all 
(2005) and Laing et al. (2006), the inner and flaring regions of the 
sources studied are defined in terms of emissivity, but the authors 
do not invoke the presence of a shock in the transition from one 
region to the next. In this paper we identify that transition as due to 
the presence of a recollimation shock. Within this frame, we com- 
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Figure 9. Logarithm of temperature (top panel) in K and logarithm of pressure (bottom panel) in Pa at the last frame of the simulation in this work. The 
position of the first recollimation shock and the terminal (or bow) shock are indicated in the pressure map. Coordinates are in parsecs. 



pare the properties of the observed discontinuity with those of its 
counterpart in the numerical simulation at z ~ 2 kpc. 

LB02b give a speed of the jet of 0.87 c in the inner region, that 
we used as the injection speed of the jet in our simulation. However, 
the flow in the simulation expands adiabatically and accelerates in 
the decreasing density ambient medium, entering the first recolli- 
mation shock with a larger speed (~ 0.98 c). The technique used in 
LB02a that allows for a fitting of the velocity of the jet is not ade- 
quately constrained in the inner region, so the initial velocity used 
here may be in error. The authors report that jet/counter-jet ratio 
is slightly smaller in the inner region than at the start of the flar- 
ing region, which could be interpreted as due to an acceleration of 



the flow in this region, as observed in our simulation. Observations 
with higher transverse resolution and sensitivity could potentially 
constrain any acceleration in the inner region. 

The position and velocity of the flow upstream/downstream of 
the standing shock depend on the acceleration of the jet in the ini- 
tial adiabatic expansion phase. Similarly, the properties of the flow 
upstream the shock determine the properties downstream the shock 
and ultimately the jet disruption. In order to study the influence of 
jet pressure and injection velocity on the evolution of the jet, we 
performed three additional simulations with the following proper- 
ties: a) Simulation 2, tij = 0.5 c, Pj /P a ,c = 3.8; b) Simulation 
3, Vj = 0.6 c, Pj/Pa,c = 7.8, and c) Simulation 4, Vj = 0.5 c, 
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Figure 10. Different profiles of variables on the jet axis at the end of the simulation. The different panels show the rest mass density (a), pressure (solid line) 
and original ambient medium pressure on the axis (dotted line; b), jet mass fraction (c), flow Lorentz factor (d), axial flow velocity (e) and temperature (f). The 
location of the first recollimation shock is indicated in the pressure plot (panel b) and the location of the disruption point is indicated in the jet mass fraction 
plot (panel c). 



Pj/Pa,c = 1. The change in the flow injection velocity and the jet 
pressure in these additional simulations produces a change in the 
power of the jet. In the case of simulation 2, the change in the in- 
jection parameters produces a reduction of the jet power of about 
an order of magnitude with respect to the original simulation. In the 
case of simulations 3 and 4, the reduction in power with respect to 



the first simulation is of a factor of 3 and 30, respectively. Finally, in 
these simulations we have used approximations to the Bessel func- 
tions that are solved in the equation of state -see the Appendix- 
dde Berredo-Peixoto et ai] |2005l :ls ervic in order to reduce 

the computational load. 

In simulation 2, the jet flow is accelerated in the pressure gra- 
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Figure 11. Profiles of radially averaged variables, weighted with the tracer and counting only those cells where the axial velocity is significantly larger than 
zero (v z > 10 -3 c). The different panels show the rest mass density (a), pressure (solid line) and original ambient medium pressure on the axis (dotted line; 
b), Mach number (c), and Lorentz factor (d). 



dient up to Vj ~ 0.82 c. The first recollimation shock occurs at 
z ~ 700 pc, closer to the source than given by LB02b, due to the 
smaller initial pressure ratio at injection and the smaller injection 
speed. After this shock, the jet is strongly mass loaded and decel- 
erated to speeds Vj < 0.5 c. The entrained jet expands and acceler- 
ates slightly to Vj — 0.6 c, but the latter expansion ends up in fur- 
ther mass loading and deceleration of the jet at z ~ 1.1 kpc. Mass 
loading and deceleration continue downstream up to z = 4.8 kpc, 
where the velocity of the jet is already subrelativistic and subsonic. 
The simulation was stopped when the head of the jet had reached 
this distance, at time t ~ 6.6 10 6 yrs. 

The jet in simulation 3 goes through the first recollimation 
shock at z ~ 800 pc. The flow enters the shock with a speed Vj ~ 
0.94 c. After the shock, the jet decelerates and is disrupted at z ~ 
1.5 kpc. As in simulation 2, mass loading continues downstream. 
The simulation was stopped at t ~ 6 10 6 yrs, when the jet head is 
at z ~ 5 kpc and the bow shock is at z ~ 6 kpc. 

Simulation 4 differs from simulation 2 in that the jet is in pres- 
sure equilibrium with the ambient medium at injection. This differ- 
ence prevents the formation of strong shocks. The jet velocity os- 
cillates close to the injection value (Vj = 0.5 c). Successive expan- 
sions and contractions, i.e. smooth recollimations, in the pressure 
gradient cause the pinching of the jet. The growth in amplitude of 
this pinching as it couples to a Kelvin-Helmholtz instability mode 



causes entrainment and the disruption of the jet. The jet is more 
stable in this case simply because the amplitude of the recollima- 
tion shocks is reduced. Nevertheless, the jet is light and slow, which 
makes it a good candidate for disruption due to the growth of in- 
stabilities jPerucho et al]|2004 120050 . At the end of the simulation 
(i ~ 7 10 6 yrs) the jet head is at z ~ 4.8 kpc. 



4 DISCUSSION 



IScheck et all d2002h performed a series of simulations of the long 
term evolution of jets with different compositions evolving in a uni- 
form density ambient medium. The initial power of those jets was 
typical for FRII sources. Compared to their simulations, the jet in 
the simulation presented here is 100 times weaker. The bow shock 
propagates at a slightly larger mean velocity in our simulation. This 
fact is due to the simulated jet being overpressured and propagat- 
ing through a decreasing density atmosphere. The m orphology of 
our jet is close to that of model LH (leptonic, hot) in lScheck et all 
(2002), although in our case, the jet is more pinched and presents 
entrainment behind the head. It is remarkable that the jet in our 
simulation is leptonic and hot and the structur e obtained is solely 
similar to the same case in IScheck etal.ld2002h . 
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4.1 Jet advance speed and source age 

Parma et ID d2002h studied a sample of FRI radio galaxies and 
estimated their spectral ages from two-frequency data. Although 
the se measurements are s ubject to uncertainties, as pointed out 
by lLaing et al] d2006h and iKatz-Stone & Rudnickl dl997n . mainly 
due to possible confusion between jet and lobe emission, we have 
used them in order to make a rough comparison of the age of the 
source in the simu lation with those of FRIs. The ages obtained by 
IParma et alj fe002l) . claimed to be lower limits, range between 10 
and 10 8 yrs, in general larger than those of FRII jets. From these es- 
timates, the authors compute advance velocities for the lobes in the 
range 1CT 3 - 1CT 2 c. VLA observations of 3C 31 (e.g., LB02a) 
show jet material at projected distances up to 150 kpc. Consider- 
ing a viewing angle of 52° (LB02a), this implies linear distances 
of about 200 kpc. From the advance velocities found in the simula- 
tion, we can put constrains on the time that the source has been con- 
tinuously active. The continuity is supported by the fact that there 
are no emission gaps in the images. At the typical advance speeds 
found in this simulation, say 5 10~ 3 c for further advance of the 
head as an upper limit, the age of the source w ould bet > 1 8 yrs , 
in agreement with the spectral ages given by IParma et alj j2002h . 
The lower limit takes into account further deceleration of the head, 
and the curved trajectory of the plasma observed in the source at 
distances larger than ~ 25 kpc. Moreover, low-resolution images 
show emission out to 300 kpc. The uncertainty in the viewing an- 
gle at large distances from the core makes difficult the use of this 
value for an estimate. Nevertheless, this fact makes of the value 
given here a strict lower limit. In the frame of intermittency models 
for the activity in AGNs, the lower limit of the age of the source 
given here puts constraints on the possible intermittency of 3C 31: 
either the periodicity is as long as 0.01 — 0.1 times the age of the 
Universe, or it has been continuously active since the onset of its 
activity. 



4.2 Early evolution and the young counterparts of FRIs 

We have shown, in Fig. [2] that the evolution of the source in 
the compact phase is divided in two stag es: the CSO-like and the 
weak CSS-like phase. Dra ke et al] d2004h have proposed that low- 
power CSS sources can be the young counterparts of large FRI 
sources. Based on spectral aging, the authors give age estimates 
of ~ 10 5 yrs for several of these CSS sources with linear sizes of a 
few kiloparsecs, and derive expansion velocities of 0.004—0.007 c. 
These velocities are in agreement with those derived from the sim- 
ulation (see Fig. [3}. Furthermore, the age of the simulated jet in the 
weak CSS-li ke phase (t ~ 10 5 — 10 6 yrs) is of the order of that 
estimated bv lDrakeetafl J2004I) for the observed low-power CSS 
sources, and its morphology is irregular, as those authors show to be 
the case of the low-power CSS jets. Fig.|2]tells us that FRIs could 
first go through a CSO stage, characterized by a regular, expand- 
ing jet, and, after developing a strong shock due to underpressur- 
ing with respect to the ambient medium, at a distance of the order 
of 1 kpc (depending on the properties of the host galaxy and the 
jet, as we have seen in Sect. [33), develop the irregular structure 
observed in the maps for the low power CSS sources. Therefore, 
powerful CSO sources could be the young FRII sources, but, from 
our results, we predict that low power CSOs could be the young 
FRI sources. 



4.3 Cocoon temperature and emission 

The high temperature of the fluid in the cocoon (Figs. l7l9l l deserves 
some discussion. Heating of the jet plasma occurs at the stand- 
ing shocks along the jet, as shown in Fig. 11 Of . The efficiency of 
the heating depends on the strength of these shocks, which ulti- 
mately depends on the jet power. In the case of the mildly rela- 
tivistic, hot jets like those considered in the present work, this last 
quantity is dominated by the internal energy density flux (or pres- 
sure). This is why simulations 2 and 4, with values of initial pres- 
sure ratio (Pj/P c ,a) of 3.8 and 1, respectively, display weaker rec- 
ollimation shocks and the temperature of jet material hardly rises 
above the injection values. In any case, the temperatures found in 
the cocoon (T ~ 10 10 K), in combination with the resulting co- 
coon densities (n e + _ ~ 10~ 3 cm -3 ), would produce a flux at 
1 MeV of vF v ~ 10~ 19 erg cm -2 s" 1 for a source located at the 
distance o f 3C 31, quite bel ow the detection limits of INTEGRAL 
and NeXT dKino et alj2007l and references therein), and the result- 
ing spectrum would also fall belo w the detection limi ts of the X-ray 
satellite XMM-Newton. Moreover. iKino et al conclude that 

the bremsstrahlung luminosity decreases with time as t , from 
the results of a cocoon expansion that follows basically the same 
time dependence as that found in this simulation. Thus, the present 
bremsstrahlung luminosity should be 10 times smaller than that 
computed here, on the basis of the calculated age of the source 
compared to the simulated time. In this respect, IKino et al] d2007l) 
show that bremsstrahlung cooling of the cocoon is only important 
in the very first stages of the evolution (t < 200 yrs). This vali- 
dates the adiabatic treatment of the problem. We want to point out 
that the bremsstrahlung emi ssion calculated fo r the shocked ambi- 
ent medium is even smaller. IZanni et al] d2003l) performed a series 
of simulations of supersonic and underdense jets in a decreasing 
pressure atmosphere and showed that jets evolve in two different 
phases regarding their high energy emission: a phase in which the 
shell formed by shocked material is highly overpressured and ra- 
diative, and a later phase in which the shock is weaker and a deficit 
of X-ray emission is expected from the lobes. The jet in our sim- 
ulation is in the former of the two phases. However, IZanni et ail 
d2003l) point out that jets with low density ratios, as that simu- 
lated here, form wide and not very dense shells from which no 
strong emission is expected, in agreement with the results given 
above. Fi nally, the recent discovery of bow sh ocks in low power 
radio jets dKraft et alj2003l : ICroston et alj2007h . moving at similar 
Mach numbers as those obtained here, and showing overpressure 
by more than an order of magnitude with respect to the ambient 
medium, gives support to our results regarding the dynamics of the 
bow shock (see next subsection). 



4.4 Bow shock 

At the end of the simulation, the head of the bow-shock has reached 
a distance of ~ 15 kpc from the injection position. It expands self- 
similarly and at basically constant rate (~ 7 10 -3 c in the axial 
direction), with a slight deceleration with time. The bow shock is 
still supersonic by the end of the simulation (M ~ 2.5), contrary to 
what is expected for an FRI jet in theoretical models, which predict 
trans-sonic speeds for shocks at s uch d istan ces. However, recen t 
X-ray observations bv lKraft et al]d2003l) and lCroston et all ( 2007) 
show the presence of bow shocks with Mach numbers between 3 
and 8 in the low power radiogalaxies Centaurus A and NGC 3801 at 
distances of a few kiloparsecs from the source. Taking into account 
that the jets in 3C 31 are among the most powerful (10 44 erg/s)in 
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FRI sources, it is plausible that the bow shock is long lived in this 
source. Furthermore, as we have shown in the previous paragraph, 
we have only simulated less than 10% of the real evolutionary time 
of the source. Thus, it is possible that, after the simulated time, 
the bow shocks may naturally decelerate to trans-sonic velocities. 
The properties and morphology of the jet would then be modified, 
accommodating the further evolution of the flow to the observed 
structure of 3C 31, where no bow-shock is observed, and the emis- 
sion appears to fade gradually with distance. 

Table [3] shows the values of number density, temperature and 
pressure at both sides of the bow shock in our simulation, com- 
pared to those given in lKraft et al] d2003l) and lCroston et alj fc007h 
for Cen A and NGC 3801, respectively. The jump in density at 
the head of the bow shock, in the axial direction, has a value of 
3.7, close to the Rankine-Hugoniot jump condition limit for strong 
shocks, similar to the case of NGC 3801. The jump in pressure in 
the simulation is smaller than in the observations of Centaurus A 
and NGC 3801, which can be understood in terms of the age of the 
jets: the jet in the simulation is older than those in Centaurus A and 
NGC 3801 and therefore the bow shock is slightly closer to equi- 
librium with the external medium. The main differences between 
the simulation and the observations arise in the temperature. In a 
direction transverse to the jet axis, the temperature of the shocked 
ambient medium is of the order of that in Centaurus A. However, 
in the head of the bow shock, the temperatures of the shocked gas 
reach values of 10 9 K, too high compared with observations. This 
difference may be due to 3C 3 1 being more powerful than the ob- 
served sources and also due to the lack of cooling mechanisms in 
the simulation. Apart from this, the results obtained from the sim- 
ulation are in fair agreement with the X-ray observations of the 
aforementioned sources. Therefore, we conclude that 3C 31 possi- 
bly went through a similar stage to that observed for Centaurus A 
and NGC 3801, and only in later stages than those simulated here, 
the bow shock decelerated to transonic speeds and disappeared. 

4.5 Jet structure and the Laing & Bridle (2002a,b) model 

The final structure of the jet is analyzed in the following para- 
graphs. We observe a fast adiabatic expansion when the jet leaves 
the galactic core, as the jet propagates through the steep density 
gradient of the galaxy. This expansion ends at z ~ 2 kpc and it 
is followed by a sudden recollimation. This recollimation gener- 
ates a new overpressuring of the jet, and thus, a second expansion. 
The smoother pressure and density gradient of the atmosphere in 
the region where the second expansion takes place, causes it to ap- 
pear smoother than the first expansion. This process ends, as in the 
previous case, in a recollimation shock, at z ~ 3 kpc. The third 
standing shock is at z ~ 4.5 kpc, after which the jet, strongly en- 
training and decelerated by mass loading, generates a wide shear 
layer. 

The structure of the jet at the end of our simulation is to be 
compared with that given in the dynamical model for the jet of 
3C 31 in LB02b. The main caveat for comparison between the sim- 
ulation and the observations and modelling in the latter paper is that 
the jet in our simulation has not reached a steady-state, but keeps 
evolving. This is clearly observed in the overall structure of the jet 
seen in Figs. |7|9| compared to the images of the jets in 3C 31. The 
main differen ce is the presence of lobes in the simulation: these are 
not observed. lLaing et al.1 feOOd) have shown that the kinematics 
of the jets in 3C296 (which are surrounded by their lobes) may be 
affected by the absence of a shear layer with the ambient medium. 

In LB02b, the inner 12 kpc of the jet are divided into three 



regions: inner, flaring and outer. The inner region (up to 1.1 kpc) 
consists of a fast moving flow v ^ 0.8 — 0.9c that enters the flar- 
ing region (1.1 — 3.5 kpc) through a discontinuity that decelerates 
the flow and increases the emissivity. In this region, the authors ob- 
serve a spread of the isophotes and a subsequent recollimation. The 
outer region (3.5 — 12 kpc) is characterized by a slow decrease in 
velocity, and continuous mass-loading. In our simulation, the dis- 
continuity between the inner and flaring region is identified with 
the recollimation shock at z ~ 2 kpc. LB02a argue that the veloc- 
ity of the jet is about 0.8 c up to about 3 kpc. In our case, however, 
after the standing shock, the velocity already drops to about 0.4 c. 
Moreover, LB02b found that the entrainment rate required to coun- 
terbalance the effects of adiabatic expansion and keep the velocity 
fairly constant at the beginning of the flaring region was consistent 
with that expected from mass injection by stellar winds. This points 
towards a fundamental difference between our simulation and their 
modelling, mainly due to the lack of mass load from stellar winds 
and the presence of a standing shock in our case. However, such a 
shock could explain discontinuities in emissivity found close to the 
start of the flaring regions in the observations and m o delling of the 
source s stu died by LB02a.b. | Canvin & Laing] J2004l) . lcanvin et all 
d2005h and lLaing et al.1 ( 2006), and the obser vation of h igh energy 
emission f rom the flaring reg ions in 3C 31 (Hardcastle et al. 2002), 
NGC 315 dWorrall et alj|200% and 3C 296 faardcastle et alj|2005l : 
lLaing et all 20061) . Onthe other hand, strong recollimation is not 
evident from the shapes of the jets in the sources studied, although 
it could be occurring in the inner parts of the jet, whereas the outer 
layers show constant or increasing radius, as shown by our sim- 
ulation, where the radius of the outer layers shows little signs of 
recollimation, contrary to the case of the core of the jet (Fig. |4p 
and Figs.[7l9t. 

In the simulation, the transition between the flaring and outer 
regions occurs at z ~ 4 kpc, after the third recollimation shock, 
and not as a slight underpressuring and recollimation of the jet that 
generates a smooth continuation between both regions, as required 
by the model of LB02b. Although there is no observational evi- 
dence for recollimation shocks in this region, these outer shocks 
in the simulation are milder than the first, so no strong emission 
from such structures would be expected. The oscillations of the jet 
pressure around pressure equilibrium captured by the simulation in 
the outer region are still strong enough to generate shocks, though 
milder than the previous ones. This difference between the simula- 
tion and the model reflects the disagreement between the assump- 
tion of pressure equilibrium at long distances to the core in LB02b 
model and the numerical results that display an overpressured co- 
coon by the end of the simulation. 

In order to illustrate the discussion in the two previous para- 
graphs, Fig.ll2lshows the jet radius versus distance. The plot clearly 
shows three regions for the jet morphology as in LB02a,b. The tran- 
sitions between the inner and flaring regions and between the flar- 
ing and outer ones are indicated in the plot. The initial expansion in 
the inner region, expansion -due to the sudden increase in pressure 
and density behind the recollimation shock- and contraction in the 
flaring region, and further expansion of the jet in the outer region, 
are observed here, as pointed out by LB02a (see their Fig. 4). The 
oscillations of the jet radius in the outer region are due to the irreg- 
ularities of the flow. Dotted lines show linear fits of the radius of 
the jet with respect to distance, that allow to give estimates on local 
jet opening angles: 13.8° for the inner region, 4.8° for the first half 
of the flaring region, and 17.5° for the outer one. These values dif- 
fer from those obtained in LB02a (approximately 8.5°, 18.5° and 
13.1° for the inner, flaring and outer regions, respectively). The 
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Table 3. Shell (shocked ambient medium in the simulation) and ISM (unperturbed ambient medium in the simulation) thermodynamic values in the simulation, 
compared to those of Centaurus A and NGC 3801. The ranges in the values given for the simulation stand for the different values of the parameters in the 
transversal and in the axial directions: the maxima (minima) in the shell correspond to the axial (transversal) direction and the maxima (minima) in the ISM 
correspond to the transversal (axial) direction. The ranges in the other columns are taken from the referenced papers. 





Simulation 


Centaurus A 1 


NGC 3801 2 


n 3 hell (™ 3 ) 


(3.3 - 5.2) I0 4 


2.0 10 4 


(2.0 - 3.0) 10 4 


niSM (m~ 3 ) 


(1.4 - 1.5) 10 4 


1.7 10 3 


4.6 10 3 


PshM (Pa) 


(0.4 - 1.8) 10~ n 


2.1 10" 11 


(4.2 - 8.9) 10" 12 


P ISM (Pa) 


1.3 10~ 12 


1.0 10~ 13 


3.8 10" 13 


Tshell (K) 


3.2 10 7 - 1.0 10 9 


3.34 10 7 


(0.8 - 1.2) 10 7 


T IS M (K) 


1.7 10 7 


3.36 10 6 


2.67 10 6 



1 Kraft et al. ( 2003). ^Croston et all 12007). 



main difference appears in the flaring region and is probably due 
to the series of recollimation shocks acting in this region in our 
case. However it should be also pointed out that we measure the ra- 
dius of the jet using velocity as a reference, whereas LB02a use the 
emitting material. We have shown in Sec. l3.1l that the jet is radially 
stratified in the simulation (Fig. [4]», so the emitting material could 
be surrounded by a non-emitting, slower wind that we consider as 
part of the jet when computing its radius (e.g., gas moving with 
velocities between 0.3 c and 0.5 c). 

The jet in the simulation accelerates to a higher velocity (see 
Fig. |10t than that given in LB02a,b for the jet in the transition from 
the inner to the flaring region. In the numerical simulation, the jet 
is injected into the grid with the speed given in LB02b for the jet 
in the inner region (0.87 c), and acceleration down the pressure 
gradient speeds up the jet to ~ 0.98 c. The differences in veloc- 
ity of the flow and position of the shocks were studied by means 
of three additional simulations with different values for velocity 
and pressure ratio with the external medium at injection: Simula- 
tion 2, with Vj = 0.5c and Pj/P a ,c = 3.8; Simulation 3, with 
Vj — 0.6 c and Pj/P a ,c = 7.8, and Simulation 4, with Vj = 0.5 c 
and Pj/P a ,c ~ 1. In simulations 2 and 3 the first recollimation 
shock occurs at z ~ 700 — 800 pc, too close to the source com- 
pared to observations, due to the smaller velocity (simulations 2 and 
3) and smaller pressure ratio with respect to the ambient medium 
at injection (simulation 2). This result points towards the values 
used in simulations 1 (main simulation) and 3 for jet overpres- 
sure (Pj/Pa,c ~ 8) and speed at injection (vj ~ 0.6 — 0.87) at 
z ~ 500 pc to be close to those in the real jet in 3C 31, as we are 
able to reproduce the transition between the inner and the flaring 
regions at about the appropriate distances, given in LB02b. Fine 
tuning of the initial jet velocity and pressure ratio with the exter- 
nal medium can certainly reproduce the exact observed position of 
the standing shock in the jet in 3C 31. As we have already pointed 
out, the significant overpressure of the jet at z ~ 500 pc is re- 
quired in order to produce a discontinuity (recollimation shock) in 
the transition from the inner to the flaring region at the observed 
position. This is confirmed by results from simulations 2 and 4, 
which have lower initial values of jet pressure: In simulation 2 the 
recollimation shock forms too close to the source compared with 
the observations, and in simulation 4 the jet shows neither signif- 
icant expansions nor strong recollimation shocks that can explain 
the increase in emission at the beginning of the flaring region. We 
have not studied in this w ork the possible i nfluen ce of a change 
in the atmosphere model of Hardcastle et al. (2002). However, any 
change in the density and pressure gradients would certainly influ- 



1200 



1000 - 




1000 2000 3000 4000 5000 6000 7000 
z (pc) 



Figure 12. Jet radius versus distance in the last frame. Jet radius is com- 
puted taking the outermost position where the axial velocity is larger than 
0.3 c. The boundaries between regions in the simulation are marked with 
dashed vertical lines; dotted lines indicate fitted parts in order to obtain 
opening angles. 

ence the position of the recollimation shock and the properties of 
the jet at this point. 



4.6 Mass entrainment 

From the panel of the tracer distribution at the end of the numeri- 
cal simulation (Fig.|7j, strong mass loading of the flow is observed 
for distances z > 4.5 kpc. At distances shorter than 4.5kpc am- 
bient material is entrained through the jet boundary. In the region 
z ~ 4.5 — 6 kpc the ambient material is entrained by the termi- 
nal shock. After this shock, the jet flow is mainly subsonic favour- 
ing the entrainment of ambient material. For comparison, in LB02b 
model, the strongest entrainment occurs at z ~ 3 — 3.5 kpc (see 
their Fig. 11) downstream the flaring point. After this local max- 
imum in the entrainment rate follows a monotonous increase be- 
yond z — 4 kpc. The authors claim that the mass entrainment is 
due to stellar mass loss near the flaring point, but due to mixing in 
a boundary layer farther away. In any case, the comparison between 
the mass entrainment rate as a function of distance to the source in 
LB02b model and that derived from our numerical simulation has 
to be considered with caution, as the simulation has not reached a 
steady state. 
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Our simulation allows to conclude that the cause of the dis- 
continuity between the inner and the flaring region and of the pro- 
cess of entrainment in the jets in 3C 31, is the generation of a 
strong standing shock due to the initial overpressure of the jet, in 
agreement with LB02b, and not to the grow th of Kelvin-Helmholtz 
instabilities to nonl inear amplitudes (e.g., iRosen & Hardee] 1 200d : 
IPerucho et all2005h . 



5 CONCLUSIONS 

We present here a hydrodynamical relativistic simulation on the 
evolution of a FRI jet, using parameters extracted from the mod- 
elling that LB02a,b made for the radio jets of 3C 31. The simu- 
lated jet is purely leptonic and propagates in a d ecreasing density 
and pr essure atmosphere with the profiles given in lHardcastle et al] 
d2002h for 3C 31. The simulation was followed up to t = 7.26 10 6 
yrs (about one tenth of the estimated lifetime of 3C 31). 

The simulation shows that the source can go through a weak 
CSS-like phase in the early evolution. The expansion speed, ages, 
sizes, and irregular morpholo gies are in agreement with those 
claimed by Dra ke et alj ( 120041) for low power CSS sources. The 
further evolution of the simulated jet into an FRI morphology sup- 
ports their claim that weak CSS sources are the progenitors of 
large scale FRI jets. Estimates of the age of the source derived 
from the advance velocities measured in the simulation lead to ages 
t > 10 8 yrs, in ag reement with recent studies of FRI radio galaxies 
(Par ma et al . 2002). The lack of emission gaps in the images points 
towards 3C 3 1 having been continuously active since the triggering 
of its activity phase, with no shorter-term periodicity. 

At the end of the simulation, the head of the bow shock has 
reached a distance of ~ 15 kpc propagating at an almost con- 
stant speed of 710~ 3 c up to t = 4.5 10 6 yrs, decelerating after- 
wards. The region encompassed by the bow shock can be divided 
in two parts: the cocoon, fed with jet material, hot and light; the 
shocked ambient medium region, cooler and denser. This struc- 
ture is reminiscent of the cavity/shell structure characteristic of 
the shocked regions surrounding powerful jets. The pressure in the 
shocked region decreases with time with a steeper slope than in the 
Begelman & Cioffil d 19891) model. A simple generalization of this 
model that takes into account the evolution of the shocked regions 
in decreasing density atmospheres has been used to explain the fast 
decrease of pressure as well as the self-similar expansion of the 
shocked region. Moreover, assuming self-similarity also for the co- 
coon evolution, our model is able to explain the constant character 
of temperature in the cocoon with time. 

At the end of the simulation, the bow shock is still slightly su- 
personic (Mach number ~ 2.5). This result is supported by recent 
observations of bow shocks with Mach num bers between 3 a nd 8 
in the low p ower radio galaxies C entaurus A dKraft et alj2003h and 
NGC 3801 dCroston et al]|2007l) . The fact that the jet power used 
in our simulation for the jets in 3C 31 is in the upper part of the 
range appropriate for FRI sources could explain the presence of a 
bow shock in the simulation at larger distances from the galaxy than 
observed in these galaxies. We show that the pressure and number 
density jumps acro ss the bow shock in the simulation are consis tent 
with those given in lKraft et all d2003l) and lCroston et alj d2007l) for 
Centaurus A and NGC 3801, respectively. From this fact, we con- 
clude that the jet in 3C 3 1 possibly went through the same stage as 
that observed for those sources, and that the deceleration and dis- 
appearance of the bow shock should occur at later times than those 
simulated here. A possible overestimate of the bow shock pressure 



and velocity could be caused by the lack of radiation cooling in the 
simulation. 

Focussing on the jet structure and dynamics, the simulation 
reveals that the jet expands in the pressure gradient and under- 
goes several subsequent recollimation and expansion processes that 
end up in the formation of a wide shear layer, mass loading and 
complete disruption of the flow. The parameters used in this sim- 
ulation succeed in explaining the general structure of the jet, as 
modelled by LB02b. However, the exact locations of transition be- 
tween model regions in LB02b are not reproduced here. The inter- 
nal shocks in the simulated jet are formed at different positions than 
given by the model. Also, the model predicts only one discontinuity 
in the transition between the inner and flaring regions, and a smooth 
transition from the flaring to the outer region, whereas we find two 
more discontinuities in the jet, the latter being the transition from 
the flaring to the outer region in the simulation. This difference may 
be caused by the assumption of pressure equilibrium at long dis- 
tances (z ~ 12 kpc) made in the modelling of LB02b. It should be 
also kept in mind that our simulation has not reached a steady state. 
This assumption is validated by the fact that the jets in 3C31 show 
no lobes, although these jets could have been overpressured as they 
went through the initial stages simulat ed here, as it happ ens with 
the jets in younger sources like 3C 296 dLaing et"al]|2006t) . The ex- 
tra discontinuities captured by the simulation are a consequence of 
the oscillations of the jet around pressure equilibrium with the ex- 
ternal medium, that turn out to be stronger than given by LB02b. 
Fine tuning of the injection velocity (t)j ~ 0.6 — 0.87 c) and a 
significant overpressure of the jet (as that used in our simulation, 
Pj/Pa,c ~ 8) are required at injection in the grid (z ~ 500 pc) 
in order to fit the first recollimation shock to the position given in 
LB02b. We conclude that a standing shock formed due to the initial 
overpressure of the jet is the cause for the discontinuity found be- 
tween the inner and the flaring regions and for the mass entrainment 
of the jets in 3C 31. 

The simulation presented here c onfirms the paradig m for the 
evolution o f FRI's dBicknelll 1 1984 iKomissarovl Il990dlbl iLaind 
1 199311 19961) . in which the adiabatic expansion of an overpressured 
jet is followed by subsequent recollimation in shocks and expan- 
sion processes. The jet is decelerated to transonic and subsonic ve- 
locities in these shocks, and a mixing layer is formed that finally 
disrupts the flow. However, recent X-ray observations and this sim- 
ulation indicate the existence of bow shocks still in kiloparsec scale 
FRI jets. Further two and three-dimensional magnetohydrodynamic 
simulations, including radiative cooling and realistic mass entrain- 
ment from stellar mass losses should be performed in order to dis- 
entangle the role of these processes in the structure and dynamics 
of the FRI jets. It is important to note that these simulations can be 
crucial for a deeper understanding not only of the evolution of jets 
in FRI radio galaxies, but also of their impact on the interstellar and 
intergalactic media, e.g., through heating. 
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APPENDIX A: THE EQUATION OF STATE 

The equati on of state of a relativistic perfect gas can be written in 
the form JSvngelll957l : lFalle & Kormssarov|[i996h : 



w — nimiG(£i), 
i=i 



(Al) 



(A2) 



where, w — ph, with p the proper rest-mass density and h the 
specific enthalpy, p is the pressure, ni is the number density of a 
given family of particles with mass mi, 



mi 



(A3) 



and 

In the latter equations fcs is the Boltzmann constant, T, the tem- 
perature and K v (£) are the modified Bessel functions: 



(A4) 



MO 



exp(— £ cosh 6) coshvO d6. 



(A5) 



The adiabatic exponent is derived from the definition of sound 



speed 
2 _ 1 (dp' 

-h\T P J m ' 

and turns out to be: 



r = 



Ef =1 ^(G'(6K! + i) : 



(A6) 



(A7) 



where G' represents the derivative of the function with respect to 
its argument, £/. 

In our case, we deal with two species of particles: leptons 
(electrons and positrons) and baryons (protons). From the leptonic 
and total proper rest-mass densities (p and pi, respectively), charge 
neutrality allows to obtain the proton rest-mass density, p p , and the 
corresponding number densities, ni and n p (and mass fractions, Xi 
and X p ). Then, 
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w = PiG{£,i) + PpG(^ P ) 



(A8) 



and 



(A9) 



Now, these two last equations together with the equations 
defining the set of conserved variables, eqs. {6]|-{9]|, form an im- 
plicit system from which the values of p, pi, p, w, T and the two 
components of flow velocity, v R and v z , can be derived. With this 
purpose, an iteration in the temperature is performed at each nu- 
merical cell in every time-step. 



Finally, 



(A10) 




